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Abstract Pulsar Wind Nebulae (PWNe) act as calorimeters for the relativistic pair 
Cn , winds emanating from within the pulsar light cylinder. Their radiative dissipation 

in various wavebands is significantly different from that of their pulsar central en- 
nj , gines: the broadband spectra of PWNe possess characteristics distinct from those 

hH ' of pulsars, thereby demanding a site of lepton acceleration remote from the pulsar 

magnetosphere. A principal candidate for this locale is the pulsar wind termination 
shock, a putatively highly-oblique, ultra-relativistic MHD discontinuity. This paper 
summarizes key characteristics of relativistic shock acceleration germane to PWNe, 
O ' using predominantly Monte Carlo simulation techniques that compare well with 

^ , semi-analytic solutions of the diffusion-convection equation. The array of potential 

5G ' spectral indices for the pair distribution function is explored, defining how these 



depend critically on the parameters of the turbulent plasma in the shock environs. 
Injection efficiencies into the acceleration process are also addressed. Informative 
constraints on the frequency of particle scattering and the level of field turbulence 
Qs ■ are identified using the multiwavelength observations of selected PWNe. These sug- 

Q^ , gest that the termination shock can be comfortably invoked as a principal injector 

^^ ■ of energetic leptons into PWNe without resorting to unrealistic properties for the 

^^ , shock layer turbulence or MHD structure. 

I> 

o 
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1 Introduction 

Pulsar wind nebulae (PWNe) have fascinated astronomers ever since the discovery 
;_j ■ of the Crab Nebula. This source provides the template for PWN studies because 

5^ I of the excellent multiwavelength spectral information (de Jager & Harding 1992; 

Atoyan & Aharonian, 1996; Abdo et al. 2010b) and stunning spatial imaging af- 
forded by radio (historic), optical (Hester et al. 1995 for Hubble) and X-ray (Weis- 
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skopf, et al. 2000 for Chandra) observations. Its unparalleled observational quality 
is driven by the exceptional powerhouse at its center, the high spin-down power 
Crab pulsar The central pulsar fuels the dissipation in its surrounding PWN (exem- 
plified in the seminal Rees & Gunn, 1974, and Kennel & Coroniti, 1984, models), 
with the nebula serving in a symbiotic relationship as the calorimeter for the pulsar 
over its entire spin-down history. Therefore, the interface between the central engine 
and the nebula must play a principal role in setting up the emission seen in PWNe. 
This boundary is the pulsar wind termination shock (PWTS), where the wind is 
abruptly slowed by the ram pressure of the circumstellar material; it forms the focus 
of this perspective on lepton acceleration in PWNe. 

This shock is a natural site for the acceleration of particles that spawn the non- 
thermal radiation in PWNe that we observe. It should possess turbulent electro- 
dynamic fields that can energize and stochastically diffuse charges extremely effi- 
ciently. Yet, the PWTS is not the only possible site for leptonic acceleration. Field 
reconnection in and near the quasi-equatorial current sheet between the pulsar light 
cylinder and the termination shock is an alternative (e.g. Lyubarsky & Kirk 2001; 
Kirk & Skjaeraasen 2003; Petri & Lyubarsky 2007). In compact regions such as 
X-points in the striped wind, magnetic reconnection can release large amounts of 
energy as field tension is converted to heat of particles. Such a prospect needs fuller 
exploration from a theoretical standpoint. Reconnection theory needs to make robust 
predictions of distributions and injection efficiencies (from thermal gas) of acceler- 
ated populations in order to connect effectively to PWNe observations. The under- 
standing of shock acceleration is more developed in this regard, and accordingly is 
the focus of this paper. We note that the solar corona may prove a powerful test- 
ing ground for honing models of reconnection in the same way that the solar wind 
has demonstrated the general viability of diffusive acceleration at non-relativistic 
shocks. It should also be remarked that all escaping pulsar wind leptons impact the 
surface of the termination shock, whereas perhaps only a minority of such thread the 
environs of the current sheet reconnection region. Notwithstanding, reconnection in 
the near wind zone may contribute significantly to the evolution of the global MHD 
structure and associated wind parameters, as well as generate some pre-acceleration, 
both of which in turn influence the cumulative contribution of the PWTS as an in- 
jector to a pulsar wind nebula over its active lifetime. 

This paper summarizes the key aspects of diffusive acceleration at relativistic 
shocks in general, and pulsar wind termination shocks in particular. As the injec- 
tor of ultra-relativistic leptons, and ions, into PWNe, this process is only indirectly 
probed by radiation observations of nebulae. The volumetric extension of PWNe en- 
compasses significant spatial stratification of both the nebular magnetic field and the 
fluid flow speed, the model template for which is the spherically symmetric Kennel 
& Coroniti (1984) contribution. Moreover, temporal evolution is significant, with 
high energy electrons cooling rapidly over the lifetime of a PWN like the Crab, driv- 
ing synchrotron "burn-off" that is probed in the X-rays (see the review of Gaensler 
& Slane, 2006, for an extensive discussion of PWN observations and guiding in- 
terpretative material). Yet multiwavelength coverage, from radio to X-ray to high 
energy gamma-rays provides substantial constraints on the PWTS acceleration pro- 
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cess. Presuming nebular fields in the range of B ^ 0.1 mG implies pair Lorentz fac- 
tors in the range fe ^ 10^ — lO'" for the Crab to enable y-ray synchrotron emission. 
Since the pulsar is unlikely to generate such energetic particles, this demands effi- 
cient acceleration at the PWTS or elsewhere. Other PWNe impose similar require- 
ments. The PWTS energy budget divides into three components: (i) thermal down- 
stream heat, (ii) turbulent fields, and (iii) non-thermal shock-accelerated leptons, 
and perhaps ions. The balance between these is not yet fully understood, though 
indications from plasma simulations are that these components are not widely dis- 
parate in their energy densities. The central acceleration issue for PWN studies is 
whether a quasi-perpendicular termination shock can generate a sufficient injection 
efficiency £inj and the right spectral index in different energy ranges. In this paper, 
it becomes evident that the index issue can be satisfied in global terms using the cur- 
rent understanding of diffusive acceleration at relativistic shocks, while more work 
is needed to address the injection issue in a satisfactory manner 



2 Lepton Acceleration at Relativistic Shocks 

To understand the nature of relativistic lepton injection into the PWN, it is insightful 
to explore the general nature of particle acceleration at relativistic shocks. The key 
characteristic that distinguishes relativistic shocks from their non-relativistic coun- 
terparts is their inherent anisotropy of the phase space distribution function /(p) at 
any position. This is due to rapid convection of particles through and downstream 
away from the shock, since particle speeds v are never much greater than the down- 
stream flow speed U2 ~ c/3 : particle distributions never realize isotropy in either 
fluid or shock rest frames. This renders analytic approaches more complicated (Pea- 
cock 1981) than in non-relativistic systems. Early analytic offerings on particle ac- 
celeration at relativistic shocks focused on solutions of the diffusion-convection dif- 
ferential equation in the test-particle approximation (e.g.. Kirk & Schneider 1987a; 
Heavens & Drury 1988; Kirk and Heavens 1989). These generally specialized to 
the limit of extremely small angle scattering (SAS, or pitch angle diffusion). In par- 
ticular, the eigenfunction solution technique of Kirk & Schneider (1987a) was later 
successfully extended by Kirk et al. (2000) to the specific case of parallel, ultrarel- 
ativistic shocks, i.e. those with upstream fluid flow Lorentz factors Fi ^ 1 in the 
shock rest frame. Kirk et al. demonstrated that as Fi — )► oo, the accelerated parti- 
cle distribution power-law index a (for dN /dp °^ p^f{p) °= p'^) asymptotically 
approached a constant, a -^ 2.23, a value realized when Fi kJ 10. This result has 
been popularly invoked in astrophysics models of various sources, but is of very 
restricted applicability, as will become evident below. While diffusion-convection 
differential equation approaches are usually restricted to SAS that would be appli- 
cable to particle transport in quasi-linear field turbulence regimes, recently they have 
been generalized by Blasi & Vietri (2005) and Morlini, Blasi & Vietri (2007) to in- 
corporate large angle deflections in MHD turbulence of larger amplitudes 8B/B . 
The operating definition of such large angle scattering (LAS) is that the particle 
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experiences momentum deflections on typical angles 0scatt ^ 1 /Ti in interactions 
with MHD turbulence in the shock environs. Clearly, for ultra-relativistic shocks, 
LAS can be realized with quite modest deflections. 

A central limitation of these analytic methods is that they are restricted to power- 
law regimes, which are only realized when there is no preferred momentum scale, 
i.e. far above the thermal injection momentum. Therefore they provide no probes of 
the injection efficiency gjnj (defined to be the fraction of particles by number resid- 
ing in the non-thermal tail of the distribution), how ei„j connects key shock envi- 
ronmental parameters, and therefore how it correlates to the non-thermal distribution 
index a . Hence the niche for Monte Carlo techniques for modeling diffusive trans- 
port in shocks. Such complementary simulation approaches have been employed 
for relativistic shocks by a number of authors, including test-particle analyses by 
Kirk & Schneider (1987b), Eflison, Jones & Reynolds (1990), and Baring (1999) 
for parallel, steady-state shocks, and extensions to include oblique magnetic fields 
by Ostrowski (1991), Ballard & Heavens (1992), Bednarz & Ostrowski (1998), El- 
lison & Double (2004), Niemiec & Ostrowski (2004), Stecker, Baring & Summer- 
lin (2007) and Baring & Summerlin (2009). The Monte Carlo method successfully 
reproduced the asymptotic Fi — !■ oo index value of a « 2.23 in work by different 
groups (Bednarz & Ostrowski 1998; Baring 1999; Achterberg, et al. 2001; Ellison & 
Double 2002). There are two main types of Monte Carlo simulation on the market: 
those that inject prescribed field turbulence to effect diffusion of charges (e.g. Os- 
trowski 1991; Bednarz & Ostrowski 1998; Niemiec & Ostrowski 2004), and those 
that describe the diffusion by phenomenological scattering parameters (e.g. Ellison, 
Jones & Reynolds 1990; Eflison & Double 2004; Baring & SummerHn 2009). It is 
this latter variety that will form the focus in this exposition, because of its ability to 
survey the parameter space of acceleration characteristics in an incisive fashion. 

Before outlining the essentials of the Monte Carlo technique used to generate 
many of the results presented here, it should be noted that there is a third popu- 
lar approach to modeling particle acceleration at relativistic shocks: full plasma or 
particle-in-cell (PIC) simulations (e.g. Hoshino, et al. 1992; Nishikawa, et al. 2005; 
Medvedev, et al. 2005; Spitkovsky 2008). PIC codes compute fields generated by 
mobile charges, and the response of the charges to the dynamic electromagnetic 
fields. Accordingly they are rich in their information on shock-layer electrodynam- 
ics and turbulence, but pay the price of intensive demands on CPUs. This presently 
limits them to exploration of thermal and suprathermal energies, so that full plasma 
simulations generally exhibit largely Maxwellian distributions (Hoshino, et al. 1992; 
Nishikawa et al. 2005; Medvedev, et al. 2005). However, we note the isolated recent 
suggestion (Spitkovsky 2008; Martins et al. 2009; Sironi & Spitkovsky 2009) of 
non-thermal tails spanning relatively limited range of energies, generated by diffu- 
sive transport in PIC simulations, with the thermal population still dominating the 
high-energy tail by number. To interface with astrophysical spectral data, a broad 
dynamic range in momenta is desirable, and this is the natural niche of Monte Carlo 
simulation techniques. 
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2.1 The Monte Carlo Method 

As informative background to the ensuing results on relativistic planar shocks, the 
structure of the simulation used to calculate diffusive acceleration is now described. 
It is a kinematic Monte Carlo technique that has been employed extensively in su- 
pernova remnant and heliospheric contexts, and is described in detail in numerous 
papers (e.g. Ellison, Jones and Reynolds, 1990, hereafter EJR90; Jones & Ellison 
1991; Ellison & Double 2004; Baring & Summerlin 2009). It is conceptually simi- 
lar to Bell's (1978) test particle approach to diffusive shock acceleration, and essen- 
tially solves a Boltzmann transport equation for arbitrary orientations of the large 
scale MHD field B. The background fields and fluid flow velocities on either side 
of the shock are uniform, and the transition at the shock is defined by the standard 
relativistic MHD Rankine-Hugoniot conservation relations (e.g. Double et al. 2004) 
that depend on both the sonic and Alfvenic Mach numbers. Particles are injected 
upstream of the shock with a Maxwell-Boltzmann distribution of finite temperature, 
and convect and gyrate in the laminar electromagnetic field, with their trajectories 
being governed by a relativistic Lorentz force equation in the frame of the shock. 
The upstream fluid frame magnetic field is inclined at an angle ©Bfi to the shock 
normal. Because the shock is moving with a velocity u(x) relative to the plasma rest 
frame, there is, in general, a u x B electric field in addition to the bulk magnetic 
field. Particle interactions with Alfven wave and other hydromagnetic turbulence is 
modeled by using a phenomenological scattering of the charges in the rest frame of 
the plasma. The scattering precipitates spatial diffusion of particles along magnetic 
field lines, and to a varying extent, across them as well. The scatterings are also 
assumed to be quasi-elastic, an idealization that is usually valid because in most as- 
trophysical systems the flow speed far exceeds the Alfven speed, and contributions 
from stochastic second-order Fermi acceleration are small. The diffusion permits a 
minority of particles to transit the shock plane numerous times, gaining energy with 
each crossing via the coherent shock drift and diffusive first-order Fermi processes. 
A continuum of scattering angles, between large-angle or small-angle cases, can 
be modeled by the simulation. In the local fluid frame, the time, 5tf , between scat- 
terings is coupled (EJR90) to the mean free path, A , and the maximum scattering 
(i.e. momentum deflection) angle, 0scatt via 5tf w X9g^^^^/{6v) for particles of 
speed V « c. Here the mean fee path is that for turning the particles around along 
field lines. Usually A is assumed to be proportional to a power of the particle mo- 
mentum p (see EJR90 and Giacalone, Burgess and Schwartz, 1992, for microphys- 
ical justifications for this choice), and for simplicity it is presumed to scale as the 
particle gyroradius, r^ , i.e. A = rjrg °^ p . Simulation results are fairly insensitive to 
this choice. Moreover, the scattering law is generally assumed to be identical in both 
the upstream and downstream fluids. Departures from this can easily be accommo- 
dated, but usually incur only a change in the spatial scales for diffusion either side 
of the shock. The parameter 77 in the model is a measure of the level of turbu- 
lence present in the system, coupling directly to the amount of cross-field diffusion, 
such that J] = 1 corresponds to the isotropic Bohm diffusion limit, where the field 
fluctuations satisfy 5B/B ^ 1 . In the quasi-linear regime, SB/B <C 1 , one expects 



6 Matthew G. Baring 

that 77 should scale inversely as the variance of the field, i.e. rj °^ (dB/B)^^ . In 
kinetic theory, rj couples the parallel (k"|| = Av/3) and perpendicular (k±) spa- 
tial diffusion coefficients via the relation Kj_/k^^ ~ 1/(1 + 77^) (Forman, Jokipii & 
Owens 1974; Ellison, Baring & Jones 1995). In parallel shocks, where the B field 
is directed along the shock normal (©Bfi = 0), 77 has only limited impact on the 
resulting energy spectrum, principally determining the frequency of scattering and 
hence the diffusive spatial scale normal to the shock. However, in oblique relativis- 
tic shocks where ©Bf 1 > , the diffusive transport of particles across the field (and 
hence through the shock) becomes critical to retention of them in the acceleration 
process. Accordingly, for such systems, the interplay between the field angle and the 
value of 77 controls the spectral index of the particle distribution (Ellison & Double 
2004; Baring 2004), a feature that is central to the interpretation of PWN spectra. 

It should be remarked that this phenomenological description of diffusion in 
Monte Carlo techniques is most appropriate at high energies (where it is more or 
less commensurate with results from Monte Carlo codes that inject prescribed tur- 
bulence), and omits the details of microphysics present in plasma simulations such 
as PIC codes. In the injection domain at slightly suprathermal energies, the influ- 
ences of complex turbulent and coherent electrodynamic effects become important, 
and will substantially modify the picture from that of pure diffusion that is presented 
here; such is the niche of PIC simulations. Note also that all subsequent simulation 
results presented here are obtained in the test particle approximation, where the 
accelerated population is not permitted to modify the overall MHD shock structure. 



2.2 Results for Relativistic Shock Acceleration 

Representative particle differential distributions dN /dp °^ p^/(p) that result from 
the simulation of diffusive acceleration at mildly-relativistic shocks are depicted 
in Figure [U (adapted from Baring 2009); the reader can survey Ellison & Double 
(2004), and Stecker, Baring and Summeriin (2007, hereafter SBS07) for Fi > 1 
simulation results that possess similar character to the parallel shock ( ©Bfi = 0° ) 
examples in the Figure. These distributions are obtained just downstream of the 
shock and are measured in the shock rest frame. They are equally applicable to 
electrons or ions, and so the mass scale is not specified; presuming that the wind 
loss from pulsars is dominated by pairs, the mass scale is nominally m^ . A striking 
feature is that the slope and shape of the non-thermal particle distribution depends on 
the nature of the scattering. The often cited asymptotic, ultrarelativistic index of a ~ 
2.23 for dN /dp °^ p^" mentioned above is realized only for parallel shocks with 
©Bfi ~ 0° in the mathematical limit of small (pitch) angle diffusion (SAS), where 
the particle momentum is stochastically deflected on arbitrarily small angular (and 
therefore temporal) scales. As mentioned above, in practice, SAS results when the 
maximum scattering angle 0scatt is inferior to the Lorentz cone angle l/F} in the 
upstream region. In such cases, particles diffuse in the region upstream of the shock 
only until their velocity's angle to the shock normal exceeds around l/Fi, after 
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Fig. 1 Particle distribution functions dN/dp from mildly-relativistic shocks (FijSi = 3, i.e. 
ji\ = Wi/c = 0.949) of upstream-to-downstream velocity compression ratio r = iiix/u2x ~ 3.24. 
Simulation results can be divided into two groups: parallel shock runs (©Bfl = 0° > upper three his- 
tograms), and oblique, superluminal shock cases (©Bfi = 20°, 40°, 60° , lower three histograms). 
Scattering off hydromagnetic turbulence was inodeled by randomly deflecting particle moinenta 
by an angle within a cone, of half-angle 6scatt . whose axis coincides with the particle momentum 
prior to scattering; the ratio of the diffusive mean free path A to the gyroradius r^ was fixed at 
T] = A/rg = 5 . The heavyweight lines (two uppermost histograms) are for the large angle scatter- 
ing cases (LAS: 1/Fi <g 9scatt < ^)- All other cases constitute pitch angle diffusion (sinall angle 
scattering: SAS) runs, when 9scatt ^ l/^i ^"d the distributions becoine independent of the choice 
of flscatt • AH distributions asymptotically approach power-laws dN/dp « p^" at high energies. 
For the two cases bracketing the results depicted, the power-laws are indicated by lightweight lines, 
with indices of CT = 1.61 (©Bfi = 0° , Sgcatt < J^) and (7 = 3.31 (©Bfi = 60° , Sgcatt < 10° ), 
respectively. Also displayed is an indication of the index required to match Fenni-LAI > 1 GeV 
observations for the Crab Nebula, assuming uncooled inverse Compton emission. 



which they are rapidly swept downstream of the shock. The Figure indicates clearly 
that when the field obliquity ©Bfi increases, so also does the index a , with values 
greater than (7^3 arising for ©Bfi ^ 50° for this mildly-relativistic scenario. This 
is a consequence of more prolific convection downstream away from the shock. 
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Figure [T] also shows results for large angle scattering scenarios (LAS, with 
4/ri ^ ©scatt ^ '!^), where the distribution is highly structured and much flatter 
on average than p^^ . The structure becomes more pronounced for large Fi (see 
Baring 2004; Eflison & Double 2004; SBS07, for details), and is kinematic in ori- 
gin, where large angle deflections lead to fractional energy gains between unity and 
r^ in successive shock crossings. Each structured bump or spectral segment corre- 
sponds to an increment of two in the number of shock transits (Baring 2004). For 
p :^ mc, they asymptotically relax to a power-law, in this case with index C7 « 1 .61 . 
Intermediate cases are also depicted in Figure [1] with Qscatt ~ 4/ri . The spectrum 
is smooth, like for the SAS case, but the index is lower than 2.23. From the plasma 
physics perspective, magnetic turbulence could easily be sufficient to effect scat- 
terings on this intermediate angular scale, a contention that becomes even more 
germane for ultrarelativistic shocks (SBS07). Note that there is a clear trend (e.g. 
see EJR90; Baring 2004; SBS07) of declining a for higher Fi , the consequence of 
an increased kinematic energy boosting in collisions with turbulence. 

The plot in Figure [T] includes an indication of the particle distribution index re- 
quired to match the Fermi observations of the Crab Nebula. The Fenni-LALT spectral 
index in the 1-20 GeV range, corresponding to a putative inverse Compton signal, 
is Uy = 1.64 (see Abdo et al. 2010a, and specifically Figure 5 therein). In the case 
where this corresponds to the in situ accelerated population (i.e. the population is 
uncooled on the relevant timescales), one finds that a ~ 2ay — I = 2.28. Such a 
scenario is depicted by the "Crab 1-20 GeV" line in Fig. [T] In contrast, if inverse 
Compton cooling is sufficiently rapid as to define the total Fermi > 1 GeV spectrum, 
then (7 = 2ay — 2 = 1 .28 . Thus, strongly-cooled IC models would suggest large an- 
gle scattering is active in the Crab Nebula termination shock, if it is superluminal. 

Now for an important definition pertaining to the following discussion. The MHD 
phase space of relativistic shocks bifurcates neatly into two regimes. In general, 
Monte Carlo simulations "operate" in a shock rest frame named the normal inci- 
dence frame (NIF), where the upstream flow is directed along the shock normal 
(usually chosen to be the x -direction, a convention adopted here). In this frame, 
the upstream magnetic field is inclined to shock normal by an angle of ©bsI ■ Due 
to relativistic aberration effects, generally ©bsI t^ ©Bfi , with equality arising only 
in truly non-relativistic shocks. For many systems, there is also a shock rest frame 
called the de Hoffman-Teller (HT) frame (identified by de Hoffman & Teller 1950), 
which is obtained by a boost miht = j3ihtc ~ Mi.v/cosQBfi along the magnetic 
field so as to bring the shock to rest. In this HT frame, there are no static electric 
fields, implying no E x B drifts parallel to the shock plane. Subluminal shocks are 
defined to be those where the HT flow speed j8iht corresponds to a physical speed, 
less than unity, i.e. the upstream field obliquity satisfies cos0Bfi < j3iv = mia/c 
When jSiHT > 1 , the de Hoffman-Teller frame does not exist, and the shock is said 
to be superluminal. This division naturally demarcates a dichotomy for the gyra- 
tional characteristics of charges orbiting in the shock layer Subluminal shocks per- 
mit many gyrational encounters of charges with the shock interface, and therefore 
also reflection of them into the upstream region. This implies efficient trapping (e.g 
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Fig. 2 Particle distribution functions dN /dp from mildly-relativistic sub-luminal shocks 
(FijjSi.v = 0.577, i.e. jSij = »i.,/c = 0.5) of upstream-to-downstream velocity compression ra- 
tio r = u\x/u^x ~ 4. Simulation results are depicted for two upstream fluid frame magnetic 
field obliquities, labelled by their corresponding de Hoffman-Teller frame upstream flow speeds 
J^iHT = i8i,Y/cos0gfx . These are in distinct groups of four: ©Bfi = 48.2° (jSijjT = 0.75, inul- 
tiplied by 10*) for the upper four histograms, and ©Bfl = 59.1° (jSiht = 0.975) for the lower 
four histograms. Scattering off hydromagnetic turbulence was modeled by randomly deflecting 
particle momenta by an angle within a cone, of half-angle 6scatt , whose axis coincides with the 
particle momentum prior to scattering; four different ratios of the diffusive mean free path X to 
the gyroradius r^ were adopted for each ©Bfl ■ All results were for small angle scattering (SAS), 
when 6scatt *S 1 /^i ^nd the distributions become independent of the choice of 9scatt • A low 
sonic Mach number Mg was chosen so as to effectively maximize the efficiency of injection from 
thermal energies. Adapted from Baring & Summerlin (2009). 



see Baring & Summerlin 2009), and effective acceleration. In contrast, for super- 
luminal shocks, in the absence of deflections of particles by magnetic turbulence, 
the convective power of the flow compels particles to rapidly escape downstream 
(e.g. Begelman & Kirk 1990), thereby suppressing acceleration. In such cases, par- 
ticles sliding along the magnetic field lines would have to move faster than the speed 
of light in order to return to the upstream side of the shock. Such dramatic losses 
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from the acceleration mechanism can only be circumvented by strong cross field 
diffusion precipitated by large amplitude field turbulence fields (e.g. Jokipii 1987; 
Elhson, Baring & Jones 1995), i.e. essentially close to the Bohm limit. 

The focus now turns to displaying the array of expectations for subluminal rela- 
tivistic shocks. Principally, we will explore how the distribution index a and injec- 
tion efficiency depend on the effective frequency A/r^ of scatterings, and the up- 
stream field obliquity ©Bfi . Representative particle (lepton or ion) differential dis- 
tributions dN /dp that result from the simulation of diffusive acceleration at mildly- 
relativistic shocks of speed j3ix = 0.5 are depicted in Figure |2] These distributions 
were generated for Qgcatt ^ 10°, i.e. in the SAS regime. Results are displayed 
for two different upstream fluid frame field obliquities, namely ©Bfi = 48.2° and 
0Bti =59.1°, with corresponding de Hoffman-Teller frame dimensionless speeds 
of jSiHT = J3ia7 cos ©Bfi = 0.75 and 0.975, respectively. The distributions clearly 
exhibit an array of indices a , including very flat power-laws, that are not monotonic 
functions of either the field obliquity ©bj i or the key diffusion parameter r\—X/rg. 
Fig. 12] also emphasizes that the normalization of the power-laws relative to the 
low momentum thermal populations (and hence the injection efficiency Einj ) is a 
strongly-declining function of X/rg. Quantitatively, Cinj drops from 0.1 —0.2 in 
the Bohm limit cases to less than 10^'* for A/r^ ~ 10^ when j3iht = 0.975 . This is 
a direct consequence of a more prolific convection of suprathermal particles down- 
stream of the shock that suppresses diffusive injection from thermal energies into 
the acceleration process. Such losses are even more pronounced when X/rg > 10"*, 
to the point that acceleration is not statistically discernible for j3iht > 0.98 runs 
with 10"*^ simulated particles. This property is salient for the pulsar wind nebula 
context discussed below. 

A parameter survey for diffusive acceleration at a typical mildly-relativistic 
shock is exhibited in Figure |3] where only the pitch angle diffusion limit was em- 
ployed. The power-law index a is plotted as a function of the de Hoffman-Teller 
frame dimensionless speed j3iht = j3ix/cos©Bfi . It is clear that there is a consid- 
erable range of indices a possible for non-thermal particles accelerated in mildly 
relativistic shocks. A feature of this plot is that the dependence of a on field obliq- 
uity is non-monotonic. When X/rg ^ 1 , the value of a at first declines as ©sti 
increases above zero, leading to very flat spectra. As jSiht approaches and eventu- 
ally exceeds unity, this trend reverses, and a then rapidly increases with increasing 
shock obliquity. This is the character of near-luminal and superluminal shocks evi- 
dent in Fig. 12] it is caused by inexorable convection of particles away downstream of 
the shock, steepening the distribution dramatically. The only way to ameliorate this 
rapid decline in the acceleration efficiency is to reduce X/rg to values below around 
10. Physically, this corresponds to increasing the hydromagnetic turbulence to high 
levels that force the particle diffusion to approach isotropy. This renders the field 
direction immaterial, and the shock behaves much like a parallel, subluminal shock 
in terms of its diffusive character Charges can then be retained near the shock for 
sufficient times to accelerate and generate suitably flat distribution functions. This 
defines a second core property illustrated in Fig. [3] a is only weakly dependent 
on ©Bfi when X/rg < 10. Observe that the indication of the particle distribution 
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gle diffusion), for mildly-relativistic shocks of upstream flow speed jSij = ml, /c = 0.5 , and an 
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Corapton emission is operable. Note that the indices for the jSifjT = 0.75 , 0.975 cases correspond 
to those of the distributions exhibited in Fig. [2] 



index corresponding to fer»z/-LAT observations in the 1-20 GeV range (for un- 
cooled inverse Compton models; same as in Fig. [TJ suggests low values of X/rg 
and proximity of the shock obliquity to the sub luminal/superluminal boundary. This 
inference will be developed further below. 

It is appropriate to identify briefly the reason why the distribution indices ap- 
proach a ^l for subluminal shocks when A/r^ ^ 1 , i.e. the field is almost laminar. 
The origin of the extremely flat distributions with cr ~ 1 is in the coherent effect of 
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shock drift acceleration at the shock discontinuity, discussed extensively in Baring 
& Summerlin (2009). This phenomenon is due to the energy gain of charges when 
they repeatedly encounter u x B electric fields (in frames other than the HT frame) 
in gyrations straddling the shock discontinuity. Such gains are experienced between 
episodic upstream excursions as charges more or less retain gyrophases that per- 
mit reflection from the shock for long periods of time. Reducing A /rg , and thereby 
introducing extremely modest amounts of cross-field diffusion, disrupts this coher- 
ence, removes particles from the shock layer, and steepens the spectrum. It is not 
clear that astrophysical relativistic shocks can contain such low levels of turbulence 
as to access this academically interesting regime of phase space. 

In concluding this overview of particle acceleration characteristics at relativistic 
shocks, it is noted that the results from these Monte Carlo simulations are in good 
agreement with those from other techniques, such as semi-analytic numerical so- 
lutions of the diffusion-convection equation, and also other Monte Carlo research 
initiatives. In particular, the artificially high choice of the compression ratio r = 4 
in Figures|2]and|3]was adopted to facilitate comparison with the semi-analytic work 
of Kirk & Heavens (1989). The reader is referred to Baring & Summerlin (2009) 
and Baring (2010) for more details on such simulation validation. 



3 The Quasi-Perpendicular Pulsar Wind Termination Shock 

The discussion now turns to lepton acceleration in pulsar wind termination shocks. 
While it is clear that their upstream flow speeds should be ultrarelativistic, it is un- 
clear how fast they are. The historical paradigm of upstream bulk Lorentz factors 
Fi ~ 10^ in the Crab Nebula has been promulgated from the seminal work of Ken- 
nel & Coroniti (1984). Pulsars can easily generate such bulk flows propagating out 
through the light cylinder, since the accelerating potentials in their gaps must en- 
ergize primary electrons to at least 7^- ~ 10^ — 10^ . This is true for both outer gap 
models (e.g. Cheng, Ho & Ruderman 1986; Romani 1996) or polar cap scenarios 
(e.g. Daugherty & Harding 1982; 1996) for the electromagnetic dissipation zone in 
gamma-ray puslars. If radiation reaction-limited curvature emission is what is prin- 
cipally responsible for the GeV emission seen in a host of Fermi-LAT pulsars (see 
the Fermi pulsar catalog compendium in Abdo et al. 2010b), then one can simply 
derive the relation '}^Xc/pc ^ 2eMAx/3 for emission turnovers Emax -^ 5 x 10"^ (in 
units of ffigC^ ) in the GeV band. Here Pc is the magnetic field curvature radius, 
which is some fraction of the light cylinder radius R\c = Pc/{2n) for pulsar period 
P seconds. Also, Ac = Ti/{mec) = 3.862 x 10^" cm is the electron Compton wave- 
length over 2k. With lO^cm < p^ < 10^ cm, it is inferred that primaries assume 
Lorentz factors 10^ ^ 70 ^ 10^ in a broad array of young to middle-aged pulsars. 

However, pair cascading is rife in both the polar cap and slot gap/outer gap 
gamma-ray pulsar pictures. Much of the pair creation (magnetic one photon or con- 
ventional two-photon) occurs outside the gaps containing accelerating potentials. 
Several generations of pair production ensue, precipitating large pair multiplicities 
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77± '^ 10 — lO"* (e.g. see Daugherty & Harding 1982 for polar cap realizations, and 
Muslimov & Harding 2003 for slot gap results). Furthermore, similar values are ob- 
tained by De Jager (2007), who used the TeV inverse Compton flux in PSR B1509- 
58 and PSR B 1823- 13 to infer the total electron deposition integrated over the ages 
of their nebulae, thereby acting as a calorimeter for their pulsar pair multiplicities 
(see also Bucciantini, Arons and Amato 2010 for generally higher estimates for 
77±). Simple energy conservation in the cascading process trades multiplicity for 
Lorentz factor, so that most of the emergent pairs propagating outwards from the 
gap region assume typical Lorentz factors of y± ^ 7o/i7± . This then defines fiducial 
bulk Lorentz factors for the pair flow escaping towards the distant termination shock, 
so that for nebular modeling purposes Fi '^ 10^ — 10"* may be more representative 
of the flow just upstream of the PWTS than the higher Kennel & Coroniti (1984) 
value. However, we note that since the wind is strongly magnetically-dominated at 
the light cylinder, mysteriously transitioning to a plasma-dominated flow at the ter- 
mination shock (the so-called infamous a problem), conversion of Poynting flux 
to bulk plasma kinetic energy is a distinct possibility for raising the value of Fi , 
perhaps taking advantage of magnetic reconnection in and near the current sheet. 

The obliquity of the PWTS is less subject to such debate. If the termination shock 
is a fairly regular spatial structure, it must be highly oblique or an essentially per- 
pendicular shock ( ©Bfi ~ 90° ) in the equatorial wind zone, and also at much higher 
pulsar latitudes. Within the light cylinder, this follows from the winding up of the 
field in a classic Parker spiral, just like the solar wind termination shock (e.g. see 
Bogovalov, 1999, for a discussion of MHD structure in oblique rotators). Only di- 
rections outside the pulsar polar regions can possess more radial fields that permit 
the shock to be merely oblique, or even quasi-parallel. The actual solid angle (cen- 
tered on the pulsar) portion of the PWTS that is quasi-perpendicular depends on 
the obliquity of the rotator, how the virtually rigid inner magnetospheric field mor- 
phology causally maps over to the field outside the light cylinder, and how the field 
geometry is modified by plasma loading. Yet it is in all probability large, regardless 
of whether Fi is as high as 10^ or as low as 10^ . From the MHD simulations of the 
Crab pulsar wind of Komissarov & Lyubarsky (2004), it is clear the the termination 
shock is non-spherical, being radially compressed in the polar zones. One can also 
entertain the possibility that the termination shock is slightly rippled, akin to what is 
an emerging paradigm for the solar wind shock based on the surprising magnetome- 
ter and energetic particle data acquired by the Voyager 1 and 11 spacecraft in the last 
few years. This can then permit localized regions of the PWTS to be subluminal or 
marginally superluminal. Or it can provide seeds for acceleration in a perpendicular 
shock zone from remote, but merely oblique shock environs. However, observa- 
tional support for any such a conjecture is a long way off since it requires angular 
resolutions exceeding that of Hubble and Chandra to probe such PWTS geometry in 
bright PWNe like the Crab (see Hester et al. 1995 for Hubble and ROSAT images) 
and MSH 15-52 (see Gaensler et al. 2002 for Chandra imaging). 

The content of the PWTS is generally presumed to be an electron-positron pair 
plasma. This derives from the leading models for dissipation in the pulsar magneto- 
sphere: pairs are rife therein due to the relative ease of leptons being stripped from 
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the neutron star surface. Thermionic emission is possible in pulsars with higher sur- 
face temperatures. Moreover, if sufficiently intense parallel electric fields persist in 
the atmosphere, space-charge limited ion acceleration can proceed (Ruderman & 
Sutherland 1975; Arons and Scharlemann 1979). Such a prospect drove ideas that 
young neutron stars (Blasi, Epstein & Olinto 2000) and magnetars (Arons 2003) 
could act as accelerating sources of ultra-high energy cosmic rays. In the context 
of PWNe, baryonic loading of the wind that impacts the termination shock is pos- 
sible, and inherently alters the character of the shock. Low energy charges then 
become subject to cross-shock potentials in the shock layer, since the inertial (i.e. 
gyrational) scales of the different species are widely disparate (e.g. see the discus- 
sion in Baring & Summerlin 2007). This can act to redistribute the thermal energy 
of the charges, possibly enhancing the injection and acceleration efficiency of lep- 
tons by tapping the inertia of the incoming thermal ions. Even if the pulsar wind 
is pair-dominated, it is still possible that the PWTS interface picks up ions from 
the proximate hydrogenic ejecta and feeds them into the acceleration process. Ob- 
servational constraints on hadronic contributions to PWN gamma-ray emission are 
substantial. For example, the multi-zone models of multiwavelength emission in the 
Crab nebula of Atoyan & Aharonian (1996) indicate that pion decay emission from 
PWTS-accelerated protons colliding with cold ambient hydrogen lies comfortably 
below the inverse Compton signal in the 100 MeV - 1 TeV band, and is only likely 
to be detectable at energies > 10 TeV. The flat spectrum and absence of any pion 
decay feature in the fer»z/-LAT spectrum of the Crab (Abdo et al. 2010a) strongly 
suggest that the pair component of the PWTS is the most relevant. Given that both 
environmental and neutron star-driven baryonic loading are uncertain, and the obser- 
vational mandate for treating hadronic emission in PWNe is limited, the discussion 
below will focus on pure lepton models for wind nebulae. 



4 Connecting to PWN Observations 

The emphasis now turns to making direct inferences on the pulsar wind termina- 
tion shock environment and its lepton acceleration characteristics using the multi- 
wavelength observations of nebular emission. This necessarily connects to the non- 
thermal power-law distribution indices a . For the best known and most intensively- 
studied case of the Crab, the radio spectral index is quite flat at a-y = 1.26 (e.g. 
Wright et al. 1979), the X-ray index is steeper at a^ '^ 2.1 (see Weisskopf, et al. 
2000; Atoyan & Aharonian 1996), the 1-20 GeV y-ray spectrum has Oy = 1.64 
(e.g. Abdo et al. 2010a) which slowly breaks to Uy ~ 2.5 above 1 TeV. The radio 
spectrum is not flat enough for synchrotron self-absorption, and no low frequency 
turnover that would be a signature of a minimum lepton Lorentz factor is observed. 
Accordingly, injection of pairs into the acceleration process must take place at ener- 
gies below around 3-10 GeV. These characteristics are more or less representative of 
other PWNe: the radio index generally lies around 1.3 (see Gaensler & Slane 2006), 
and is flatter than the X-ray and TeV gamma-ray spectra, a nice synopsis of which 
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is provided in the recent review of Kargaltsev & Pavlov (2010). Even if effective 
radiative cooling is invoked at the maximum pair energies generating X-ray syn- 
chrotron emission, it becomes evident from these properties that the pair injection 
spectrum is convex, ranging from a ~^ 1.5 below around 30 GeV to C7 ^ 2.3 well 
above 1 TeV (e.g. see Bucciantini, Arons & Amato 2010). As will become evident 
shortly, this is a significant constraint on diffusive acceleration at the PWTS. 

We now identify the shock conditions required to generate these detected spec- 
tral indices, by considering shocks of higher speeds than in the previous Section. 
Representative spectral index results from Monte Carlo simulation runs are exhib- 
ited in Figure |4] (from Summerlin & Baring, in preparation), the j8i, = 0.71 portion 
of which mirrors those presented in Figure |3] When the Bohm limit of X/rg = 1 
is realized, the non-thermal distribution index is approximately independent of the 
field obliquity. When the shock is superluminal, the index C7 is a rapidly increasing 
function of ©Bfi . In subluminal regimes due to the powerful convective infleunces, 
when the field is laminar and A/r^ ^ 1 , very flat spectra can be realized because 
particles can be trapped in the shock layer and shock drift acceleration is very effec- 
tive. Note also, that inefficient injection from thermal energies is then operative, as is 
exhibited in Fig.|2] The /3i.t = 0.95 indices are those taken from Fig.[T]and indicate 
a moderate increase with obliquity in the superluminal regime. Such an increase is 
tempered relative to the j3i, = 0.71 , X/r^ = 10 situation largely because X/rg and 
the compression ratio are higher. The superluminal and ultra-relativistic Fi w 10, 
©Bfi = 60° results are possibly the most representative of the PWTS. They illus- 
trate a significant sensitivity of a to X/rg, yet the indices are lower than those for 
the shocks of lower speeds j3i;f This is caused by the increased kinematic energy 
gains in shock crossings for high Fi for quasi-elastic interactions between charges 
and MHD turbulence in the shock layer. 

The spectral indices observed for the Crab nebula in radio. X-ray and gamma- 
ray wavebands, as marked in Figure |4] offer clear constraints on the shock environ- 
ment, if diffusive acceleration at the PWTS is the operable injection in PWNe. These 
can be taken to be more or less representative of the broader population of PWNe, 
though variations exist in the observational database. The flat radio (synchotron) 
spectra demand that the turbulence generate large mean free paths along B if the 
PWTS is subluminal. Shocks by their nature generate turbulence at levels that make 
this scenario unlikely (see the discussion in Baring & Summerlin 2009), disrupting 
the coherence that permits shock drift acceleration to operate prolifically. Given that 
the PWTS is very probably superluminal over most of its surface, the small angle 
scattering regime cannot supply flat enough acceleration distributions. Large angle 
scattering can though, as is evident in Fig. [T] and in Stecker, Baring & Summerlin 
(2007). This is not an unduly restrictive demand in ultra-relativistic shocks, since 
LAS is delineated by deflections 0scatt <^ l/^i ' ^i^^ i'^ is c^sy to envisage that MHD 
turbulence in such shocks can spawn scattering angles of the order of a degree or so. 
It is this scenario that is the one most probably pertinent to the 1-30 GeV leptons. 

The inverse Compton gamma-ray signal measured by the Ferm/-LAT is probing 
leptons of energies in the TeV range. The spectroscopic demands are now different: 
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Shock Acceleration: Small Angle Scattering 
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Fig. 4 Power-law indices a for simulation runs in the limit of small angle scattering, for rela- 
tivistic shocks of three different speeds. The indices are displayed as functions of the fluid frame 
field obliquity 0Bfl (contrasting Fig.|3). Simulation data for the points connected by dotted lines 
were for an of upstream flow speed jSij = mi.v/c = 0.71 , and an MHD velocity compression ratio 
r = 3.02; for these runs, obliquities 0Bfl > ^5° constitute superluminal shocks. These index re- 
sults were obtained for different diffusive mean free paths A parallel to the mean field direction, 
namely A/r^ = 1 (squares), A/r^ = 10 (triangles), A/r^ — 10^ (pentagons), and A/j-j, = IQ-' (tri- 
angles), as labelled. Data for the higher shock speed (FijSi., = 3 =^ j8i., = ml,/c = 0.949) spectra 
displayed in Fig[T](SAS only) are exhibited as circular points with dots centered therein. These are 
mostly superluminal and corresponded to A/r^ = 5. The final subset of datapoints are the three 
filled squares grouped at ©Bf i = 60° for runs with Fi j3i., = 10 ( jSi^^ « 0.995 ) with r = 3.02 (Sum- 
merlin & Baring, in preparation). These were obtained for A/r^ — 1,3,6 ranging from the bottom 
to the top. As with previous Figures, short heavyweight lines are used to indicate the approxi- 
mate spectral index a that is appropriate to match Crab Nebula spectra in different wavebands: 
radio (presumed to be uncooled synchrotron), 10-100 keV X-ray (cooled synchrotron emission) 
and Fermi-LAI > 1 GeV observations (uncooled inverse Corapton); see the text for a discussion. 



the 1-20 GeV index can be supplied by either subluminal or superluminal shocks 
with SAS operating, provided that the turbulence is not far from the Bohm limit. 
LAS is also possible, but would require highly superluminal conditions to effect 
the requisite balance between large kinematic gains in shock-layer scatterings and 



Lepton Acceleration in Pulsar Wind Nebulae 17 

rapid convective losses downstream. The X-ray spectrum samples the highest en- 
ergy electrons, in the super TeV range, that are subject to strong cooling (burn-off) 
over the nebular lifetime. Allowing for this modification, the inferences for the in- 
jected lepton spectrum at the PWTS are similar to those from the gamma-ray data. 
A broadband picture emerges that is highlighted in Atoyan & Aharonian (1996): 
the electron spectrum is convex (i.e. steepening) in the sense that a is an increasing 
function of energy. This is not difficult to accommodate using the results from shock 
acceleration theory presented here, being modeled by a modest transition from LAS 
at low energies to Bohm-domain SAS at the highest pair energies. The portions of 
the PWTS driving this energization can be either superluminal or marginally sublu- 
minal. It is not hard to envisage turbulence that is slightly stronger at smaller scales 
than larger ones that might precipitate this LAS — > SAS evolution with energy or 
Larmor radius in gyroresonant interactions. Yet in the near-term future, neither can 
observations resolve the angular scales to demonstrate such, nor can plasma simu- 
lations probe the wide dynamic ranges in lengthscales to validate such a scenario. 

The simulation results presented here are for species of a single mass, obviously 
applying to pure pair shocks. It is natural to ask whether they might differ if the 
abundance of ions is significant. The answer must be deferred to future explorations 
of diffusive acceleration in hydrogenic plasma shocks. Yet it is expected that the 
the index results should be the same unless the turbulence generation is different 
when massive species are present. At energies below 1 GeV, the gyrational scales 
of protons and electrons of a given energy differ because the protons are at most 
only mildly-relativistic. This must lead to significantly different gyroresonant in- 
teractions for e^ and p. Furthermore, at these energies, charge separation cross 
shock potentials are anticipated to play a profound role in energy exchange between 
the two species (e.g. Baring & Summerlin 2007). These two contributions should 
provide substantial differences in injection efficiency between pair shocks and hy- 
drogenic or electron-ion ones. This injection issue is clearly salient for the overall 
prediction of fluxes in different bands for PWNe. However, at energies well above 
10 GeV, the gyro-scale of a charge of a given energy is independent of its mass, so 
that to leading order, turbulence generation and diffusion characteristics should be 
similar in this domain for relativistic e^ and e — p shocks. 

Another question is whether or not the well-known non-linear spectral concavity 
encountered in non-relativistic shocks that efficiently accelerate charges (see Jones 
& Ellison 1991; Ellison & Double 2002; Baring 2004, and references therein) might 
compete with and preclude the spectral convexity that is demanded by the multi- 
wavelength observations. Such non-linear enhancements of high energy particles 
arise for distributions that have indices (7^2, where these particles supply a size- 
able portion of the total energy flux through the shock, and thereby modify the global 
MHD shock structure. While the radio observations in PWNe access this domain, 
the distribution convexity demanded by the gamma-ray and X-ray data must mute 
possible non-linear modifications, so that they should play a more minor role than 
in the non-relativistic shocks that illuminate supernova remnant outer shells. 
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5 Conclusions 

This paper has outlined the key features of relativistic shock acceleration that pertain 
to lepton injection at termination shocks into pulsar wind nebulae. This shock is the 
most popular site for such injection, because (i) all pairs emanating from the pulsar 
that travel to the nebula must transit through this interface, (ii) it should be turbulent 
and therefore an efficient injector/accelerator, and (iii) the main characteristics of 
diffusive acceleration theory at shocks are fairly well understood. While it is quite 
possible that pre-acceleration can arise in magnetic reconnection zones between the 
pulsar light cylinder and the PWTS, such seed particles can be further energized at 
the shock to the point of masking the signatures of pre-acceleration. The histori- 
cal models that developed the paradigm of the PWTS as an injector predate refined 
studies of relativistic shock acceleration over the last decade. As is evident here, 
these more recent studies support such a paradigm in being able to generate the req- 
uisite distribution indices to match the multiwavelength PWN observations without 
appealing to unlikely situations concerning turbulence in the shock layer The key 
issue that remains unresolved by theory is how efficient injection arises from ther- 
mal energies in the PWTS. Does it occur for pure pair shocks, or is some baryonic 
loading necessary to precipitate prolific energization? Or, is a pre-acceleration seed 
required to set the diffusive processes at the shock active all the way to super- TeV 
energies? Addressing such questions will require more advanced simulations and 
theoretical analyses. The answer will illuminate the overall particle budget in pulsar 
wind nebulae, balancing radiation luminosity, non-thermal particle energetics, and 
the wind power from the pulsars that drives these fascinating systems. 
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